Sets
         a2              point evaluating integrals      /1*%v%/
*Modify*
         j               unknown parameters              /j1*j3/
*End*
;
alias(j,i,jj,ii);
alias(a2,b2);
Parameters
         meanPhase1
         theta2                          unknowm parameter evaluation
         log_joint2                       log of joint distribution
         weight_1_2
;
Scalar
         pi                              pi number               /3.14159265/
;
variables
         mean_posterior
         U_matrix                        orthonormal eigenvectors of cov
         zp_sq                           standarized value for the posterior
         log_q2
         Term1_2
;
positive variables
                  precision_posterior;
Equations
         Standarized_values1_2
         orthonormal
         log_q1_eq2                       posterior factor for 1st parameter
         Term1_Eq2
;
Standarized_values1_2(j,a2)$(ord(j) le card(j))..
         2*zp_sq(j,a2)=e=
            power(sum(i,
            U_matrix(j,i)*(theta2(i,a2)-mean_posterior(i))),2)*precision_posterior(j)
            ;
orthonormal(i,j)$(ord(i) ge ord(j))..
         sum(ii, U_matrix(i,ii)*U_matrix(j,ii))=e=0+1$(ord(i) eq ord(j));
log_q1_eq2(j,a2)..
         log_q2(j,a2)=e=0.5*log(precision_posterior(j)/(2*pi))-(zp_sq(j,a2));
Term1_Eq2..
   Term1_2=e=
         sum(a2,
                 weight_1_2(a2)*
                 (sum(j, log_q2(j,a2))
                 *exp(log_joint2(a2))));
model Var_EP /Standarized_values1_2,orthonormal,log_q1_eq2,Term1_Eq2/;
*model Var_EP /all/;
Var_EP.optfile=1;
parameter cov, CPUsec2;
FILE phase2_mean /phase2_mean.csv/;
phase2_mean.PC=5;
phase2_mean.nd=7;
phase2_mean.pw=2000;
FILE phase2_cov /phase2_cov.csv/;
phase2_cov.PC=5;
phase2_cov.nd=7;
phase2_cov.pw=2000;
FILE phase2_U /phase2_U.csv/;
phase2_U.PC=5;
phase2_U.nd=7;
phase2_U.pw=2000;
Table data_preprocessing(a2,a2)
$ondelim
$include eval_points.csv
$offdelim
;
Table data_mean(j,a2)
$ondelim
$include meanPhase1.csv
$offdelim
;
Table aux_U_matrix(a2,a2)
$ondelim
$include auxU_matrix.csv
$offdelim
;
meanPhase1(j)=data_mean(j,'1');
theta2(j,a2)=sum(b2$(ord(j) eq ord(b2)), data_preprocessing(a2,b2));
log_joint2(a2)=sum(b2$(ord(b2) eq (card(j)+1)), data_preprocessing(a2,b2));
weight_1_2(a2)=sum(b2$(ord(b2) eq (card(j)+2)), data_preprocessing(a2,b2));
mean_posterior.fx(j)=meanPhase1(j);
precision_posterior.l(j)=1;
loop((i,j),
         loop(a2$(ord(a2) eq ord(i)),
                 loop(b2$(ord(b2) eq ord(j)),
                         U_matrix.l(i,j)=aux_U_matrix(a2,b2);
                 );
         );
);
zp_sq.l(j,a2)=
            power(sum(i,
            U_matrix.l(j,i)*(theta2(i,a2)-mean_posterior.l(i))),2)*precision_posterior.l(j)/2
            ;
log_q2.l(j,a2)=0.5*log(precision_posterior.l(j)/(2*pi))-(zp_sq.l(j,a2));
Term1_2.l=sum((a2),
                 weight_1_2(a2)*
                 ((sum(j, log_q2.l(j,a2)))
                 *exp(log_joint2(a2)))) ;
option NLP=conopt;
option iterlim=1e5;
option domlim=1000;
option solprint=off;
option limcol=10;
option limrow=10;
option decimals=6;
solve Var_EP max Term1_2 using NLP;
CPUsec2=Var_EP.resusd;
cov(i,j)=sum(ii, ((1/precision_posterior.l(ii))*U_matrix.l(ii,i)*U_matrix.l(ii,j)));
put phase2_mean
put CPUsec2;
loop(j,
         put mean_posterior.l(j);
);
putclose;
put phase2_cov
loop((i,j),
         put cov(i,j);
);
putclose;
put phase2_U
loop((ii,j),
         put U_matrix.l(ii,j);
);
putclose





